function [pval2s, sRet, tRet] = f_avgretNW(exRet)

T     = sum(isfinite(exRet),1);

nlag = round((4*(size(exRet, 1)/100))^(2/9));
for i = 1:size(exRet, 2)
YY    = exRet(:, i);  
iota  = ones(size(YY));
output  = nwest(YY, iota, nlag);
mRet(1, i) = output.beta;
tRet(1, i) = output.tstat;
sRet(1, i) = mRet(1, i)./tRet(1, i);
end

[pval1s, pval2s] = f_pval(tRet,T);

function [pval1s, pval2s] = f_pval(tRet,T)

nvar   = length(tRet);
pval1s = nan(1,nvar);
pval2s = nan(1,nvar);

if length(T)<nvar
    T = repmat(T,1,nvar);
end

for i=1:nvar
% ONE-SIDED
if tRet(i)<0       % [H0: mu=0 H1: mu < 0]
    pval1s(1,i) = tcdf(tRet(i),T(i)-1);
elseif tRet(i)>0   % [H0: mu=0 H1: mu > 0]
    pval1s(1,i) = 1-tcdf(tRet(i),T(i)-1);
end
% TWO-SIDED
pval2s(1,i) = 2*((tcdf(-abs(tRet(i)),T(i)-1))); % [H0: mu=0 H1: mu \= 0]
end